On the dangers of Boolean networks: 
Activity dependent criticality and threshold networks not faithful to biology 
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Spin models of neural networks and genetic networks are considered elegant as they are accessible 
to statistical mechanics tools for spin glasses and magnetic systems. However, the conventional 
choice of variables in spin systems may cause problems in some models when parameter choices 
are unrealistic from the biological perspective. Obviously, this may limit the role of a model as a 
template model for biological systems. Perhaps less obviously, also ensembles of random networks 
are affected and may exhibit different critical properties. We here consider a prototypical network 
model that is biologically plausible in its local mechanisms. We study a discrete dynamical network 
with two characteristic properties: Nodes with binary states and 1, and a modified threshold 
function with 0o(O) = 0. We explore the critical properties of random networks of such nodes and 
find a critical connectivity K c = 2.0 with activity vanishing at the critical point. 

PACS numbers: 64.60. De, 87.18.-h, 05.50.+q, 89.75.-k 
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We currently experience a revived interest in dynami- 
cal networks of nodes with binary states, driven by two 
active fields of research: modeling of molecular informa- 
tion processing networks (as, e.g., genetic networks or 
protein networks) as well as modeling of adaptive 
networks Q. These network models with binary states 
are reminiscent of artificial neural networks as studied in 
the statistical mechanics community about two decades 
ago. 

An early motivation of networks with binary node 
states <7j 6 {0, 1} was given by McCulloch and Pitts in 
1943 Q as a model for neural information processing. 
A model for associative memory constructed from such 
nodes by Hopfield in 1982 Q attracted considerable in- 
terest among physicists as it is conveniently accessible to 
equilibrium statistical mechanics methods 0^3- A sim- 
ple redefinition of weights and thresholds maps the model 
onto the mathematical representation of a spin glass with 
states <7j G { — 1,1}, which has become the usual form of 
the Hopfield model in the physics literature. The corre- 
sponding redefinition of weights and thresholds does not 
affect the functioning of the model, as its mechanism of 
an associative memory works on the redefined weights as 
well. 

In some circumstances, however, when faithful repre- 
sentation of certain biological details is important, the 
exact definition matters. In the spin version of a neural 
network model, for example, a node with negative spin 
state <jj = —1 will transmit non-zero signals through 
its outgoing weights cy , despite representing an inactive 
(!) biological node. In the model, such signals arrive 
at target nodes i, e.g., as a sum of incoming signals 
hi = Y^=i c ij tJ j- However, biological nodes, as genes 
or neurons, usually do not transmit signals when inac- 
tive. In biochemical network models each node repre- 
sents whether a specific chemical component is present 
(a = 1) or absent (a = 0). Thus the network itself is 
mostly in a state of being partially absent as, e.g., in a 



protein network where for every absent protein all of its 
outgoing links are absent as well 25[. In the spin state 
convention, this fact is not faithfully represented. 

Another example for an inaccurate detail is the com- 
mon practice to use the standard convention of the Heavi- 
side step function as an activation function in discrete dy- 
namical networks (or the sign function in the spin model 
context). The convention 0(0) = 1 is not a careful rep- 
resentation of biological circumstances. Both, for genes 
and neurons, a silent input frequently maps to a silent 
output. Therefore, we use a redefined threshold function 
defined as 



e (x) 



1, x > 
0, x < 0. 
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When studying statistical properties of ensembles of 
threshold networks with random links, these details have 
a considerable influence on the network's dynamics and 
critical properties. When simulating ensembles via net- 
works of spins <Xj G { — 1,1}, care should be taken to prop- 
erly renormalize weights and activation thresholds to en- 
sure faithful implementation of the original model with 
states <Ji G {0, 1}. However, this is frequently omitted, 
resulting in the statistics of a system of limited biological 
plausibility 

Another example where normalization and the defini- 
tion of the nodes' thresholds matters are adaptive net- 
works, currently discussed in the context of neural net- 
works [3, [IsT - flq ] . When defining local adaptive mecha- 
nisms, it is particularly important to base it on biolog- 
ically plausible definitions of nodes and circuits. While 
these mechanisms work also for spin type networks [Hj], 
such an implementation is not realizable in a biological 
context, as it would require signals over links which are 
in fact silent, due to the inactivity of their source nodes. 
An adaptive algorithm based on such correlations of non- 
activity is therefore not plausible. 

In this paper, we first define a binary threshold network 
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that does not include explicitly forbidden states in the 
context of biological examples. Then we study its critical 
properties which we find to be distinctly different from 
those of random Boolean networks [17H21] and random 
threshold networks [1, 22, 23 1. In particular, activity of 
the network now influences criticality in a non-trivial way, 
as recently observed for random threshold networks with 
bistable nodes |24| . 

Let us consider randomly wired threshold networks of 
N nodes <n € {0, 1}. At each discrete time step, all nodes 
are updated in parallel according to 



ai(t+l) = @ (fi(t)) 



using the input function 
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In particular we choose Oo(0) :— for plausibility reasons 
(zero input signal will produce zero output). While the 
weights take discrete values = ±1 with equal proba- 
bility for connected nodes, we select the thresholds 6*^ = 
for the following discussion. For any node i, the number 
of incoming links ^ is called the in-degree fc^ of that 
specific node. K denotes the average connectivity of the 
whole network. With randomly placed links, the prob- 
ability for each node to actually have ki = k incoming 
links follows a Poissonian distribution: 



K k 

p{h = k) = — ■ exp(-K), 



(4) 



To analytically derive the critical connectivity of this 
type of network model, we first study damage spreading 
on a local basis and calculate the probability p s (k) for 
a single node to propagate a small perturbation, i.e. to 
change its output from to 1 or vice versa after chang- 
ing a single input state. The calculation can be done 
closely following the derivation for spin-type threshold 
networks in ref. 23[, but one has to account for the 



possible occurrence of '0' input signals also via non-zero 
links. The combinatorial approach yields a result that di- 
rectly corresponds to the spin-type network calculation 
via Ps° ol (fc) = pf ln (2k). However, this approach does 
not hold true for our Boolean model in combination with 
the defined Theta function 9o(0) := as it assumes a 
statistically equal distribution of all possible input con- 
figurations for a single node. In the Boolean model, this 
would involve an average node activity of b = 0.5 over 
the whole network. Instead we find (Fig. [T]) that the av- 
erage activity on the network is significantly below 0.5. 
At K — 4 (which will turn out to be already far in the 
supercritical regime), less than 30 percent of all nodes 
are active on average. Around K w 2 (where we usu- 
ally expect the critical connectivity for such networks), 
the average activity is in fact below 10 percent. Thus, 




FIG. 1: Average node activity b as function of connectivity K 
measured on attractors of 10000 sample networks each, 200 
nodes. 



random input configurations will more likely consist of 
a higher number of '0' signal contributions than of ±1 
inputs. 

Therefore, when counting input configurations for the 
combinatorial derivation of p s (k), we need to weight 
all relevant configurations according to their realization 
probability as given by the average activity b. For the 
first k = 1,2,3... this yields 
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which generalizes to 



Ps{k) = - 



EM)* 



k - 1 



(5) 



using Xt = AVi ■ (f^f) and X 1 = \. 

As the in-degree k is not equal for all nodes, the expec- 
tation value of p s (K) is essential to determine the critical 
connectivity of the whole network. This will yield the 
average probability for damage spreading for a certain 
average connectivity K . 

Consider a network of size N — 1. For large N, the 
problem studied in the above section is equivalent to con- 
necting a new node j with state <r,j to an arbitrary node 
i, increasing ki to fcj + 1. If now Uj is changed, this 
will result in a state change of node i with a probability 
p s {k + 1). As mentioned above, the link distribution fol- 
lows a Poissonian; thus we can calculate the expectation 
value (p s }(K) for the thermodynamical limit N — > oo 



( Ps )(K) = exp(-K) ]T \vs{k + 1 



fc! 



(G) 



3 




n r 

-~"w - " 1/ " — — " " lit ' 

)K ,X ,."Hf 
[p * X' ,f*' 

[] m t 

[] * X if' 
1.5 e <*K 

f b 00 (K) 

[sgk f K c (b) 
1 &*- 



_i_ 



_i_ 



_i_ 



0.04 0.08 0.12 0.16 

b 

FIG. 2: Average activity b on attractors of different network 
sizes (right to left: N = 50, 200, 800, ensemble averages were 
taken over 10000 networks each). Squares indicate activity 
on infinite system determined by finite size scaling, which is 
in good agreement with the analytic result (solid line). The 
dashed line shows the analytic result for K c (b) from eq. (0. 
The intersections represent the value of K c for the given net- 
work size. 
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TABLE I: Critical Connectivity K c for different sizes as de- 
termined from curve intersections in Figure [2] 



which can be computed numerically using p s (k) as given 
in eq. (0). 

The actual calculation of the critical connectivity K c 
for the network model presented above is now done by 
applying the annealed approximation as introduced by 
Derrida and Pomeau [18| . The application of this method 
on threshold networks is discussed in detail in [23J , which 
can be directly transferred on the network model dis- 
cussed in the present work. The critical connectivity can 
thus be obtained by solving 



(p s )(K c )-K c 



(7) 



However, K c now depends on the average network activ- 
ity, which in turn is a function of the average connectivity 
K itself as shown in Fig. Q] From the combined plot in 
Fig. [2] we find that both curves intersect at a point where 
the network dynamics - due to the current connectivity 
K - exhibit an average activity which in turn yields a 
critical connectivity K c that exactly matches the given 
connectivity. This intersection thus corresponds to the 
critical connectivity of the present network model. 

However, the average activity still varies with different 
network sizes, which is obvious from Figure[2] Therefore, 
also the critical connectivity is a function of N. Table U 
lists results for different values of N. For an analytic ap- 
proach to the infinite size limit, we can now calculate the 



probability for a node at given in-degree k and average 
network activity b t at time t to exhibit output state 1. 
This probability equals the average activity for the next 
time step b t +\. By examining all relevant input configu- 
rations, we find that for given constant k this generalizes 
to 
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Again, we have to account for the Poissonian distribution 
of links in our network model, so the average evolution 
of network activity is obtained by 



(b t+1 )(K) =exp(-K)J2^rb t+ i(k). 
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(9) 



It is now possible to distinguish between the different dy- 
namical regimes by solving (&t+i) — bt{K) for the critical 
line. The solid line in Figure [2] depicts the evolved ac- 
tivity in the long time limit. We find that for infinite 
system size, the critical connectivity is at 

K C (N -> oo) = 2.000 ± 0.001 

while up to this value all network activity vanishes in 
the long time limit (b^c, = 0). For any average connec- 
tivity K > 2, a certain fraction of nodes remains active. 
In finite size systems, both network activity evolution 
and damage propagation probabilities are subject to fi- 
nite size effects, thus increasing K c to a higher value. 

As a numerical verification, we can also derive the crit- 
ical connectivity for infinite system size K C (N — > oo) 
using finite size scaling of the above simulation results 
(Table UJ). The optimum fit is shown in Figure [3] and 
yields 

K C (N ->• oo) = 2.00 ±0.01 

which perfectly supports the analytical calculation. The 
same consideration is also possible to obtain the average 
activity b(N) for increasing network size. This can be 
done both at constant values of K as well as along the 
critical line in Figure [2] using the values of K C (N) from 
Table |TJ For the critical line, we indeed find vanishing 
network activity (inset of Figure [3]) : 

6 c (A^oo) =0.00 ±0.01. 

This does also hold true for any connectivity below the 
critical line. For supercritical networks at K > K c , the 
numerical simulation yields a non-zero fraction of nodes 
which remains active on average in good agreement with 
the analytic computation (see squares in Figure [2]). 

In additional numerical simulations using the stan- 
dard step function, we obtain a critical connectivity of 
K c « 3.7, which is analytically supported by a com- 
binatorial calculation following ref. |23| where we find 
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FIG. 3: Finite size scaling of K C (N), optimum fit shown here 
for K C (N — > oo) = 2.00. Inset: Finite size scaling of b(N) 
along the critical line, the optimum fit is obtained for vanish- 
ing activity b c (N — > oo) = 0. 
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FIG. 4: Average attractor length at different network sizes. 
Ensemble averages were taken over 10000 networks each at 
(A) K = 2.4, (B) K = 2.0, (C) K = 1.6. Inset figure shows 
the scaling behavior of the corresponding transient lengths. 



Ps° ol (fc) = p s s pm (2k). The networks exhibit significantly 
higher average activity, while most of the active nodes are 
frozen in the active state. On a side note, if we chose to 
calculate the critical connectivity based on an assumed 
average activity of 0.5, the activity-dependent calcula- 
tion via ([5]) and ([7]) would effectively reproduce the same 
result (K c « 3.7) as the original combinational approach 
from [23l ] would yield for the new Boolean model. As 
the assumption does not hold true here, this result can 
only be viewed as an additional plausibility check for the 
correspondence between both approaches. In passing we 
note that also for the spin type threshold network stud- 
ied in [23] we find activity levels different from 0.5, such 
that conditions for a valid combinatorial approach might 
not be met in that case, as well. 

Finally, let us have a closer look on the average length 
of attractor cycles and transients. As shown in Fig.SJ the 



behavior is strongly dependent of the dynamical regime 
of the network. As expected and in accordance with 
early works on random threshold networks [H as well 
as random Boolean networks [20j j . we find an exponen- 
tial increase of the average attractor lengths with net- 
work size N in the chaotic regime (K > K c ), whereas 
we can observe a sub-linear increase in the frozen phase 
(K < K c ). We find similar behavior for the scaling of 
transient lengths (inset of Figured]). 

In summary we studied threshold networks with 
Boolean node states that are biologically more plausible 
than current Boolean and threshold networks and which 
are simpler than the recentl y in troduced networks with 
bistable threshold nodes 2(| 27 1 . A major observation is 
that activity of the nodes depends on connectivity which 
also renders critical properties of the networks activity- 
dependent, as found earlier for random threshold net- 
works with bistable nodes [3]. We extend the annealed 
approximation to correct for these effects and find con- 
nectivity K c = 2.0 and vanishing activity at the critical 
point in the thermodynamic limit. 

Going beyond the statistics of random network ensem- 
bles, also real biological circuits can be implemented with 
great ease using the threshold networks defined here. We 
successfully reproduced the dynamical trajectory of the 
budding yeast cell cycle network as implemented with 
bistable threshold functions in (2fJ, as well as for the 
corresponding network in fission yeast 27 1. [28[ 



To conclude, let us remind ourselves of the original idea 
of using random Boolean networks for characterizing typ- 
ical properties of biological networks [TTj. In 1969, using 
random Boolean networks as a null model for genetic 
networks was a logical approach, given our complete ig- 
norance of the circuitry of genetic networks at that time. 
Thus the best guess was to treat all possible Boolean rules 
as equally probable. Today, however, we have much more 
details knowledge about certain properties of genetic net- 
works and, therefore, about more realistic ensembles of 
random networks. A biologically motivated and carefully 
defined threshold network, as attempted in this paper, 
may provide a more suited null model for the particular 
properties of biological networks than random Boolean 
networks with equally distributed Boolean functions. 
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